#!/usr/bin/env python

import os, sys
import numpy as np
from astropy.table import Table
from scipy import odr
import matplotlib.pylab as plt

ccd_idx = sys.argv[1]
ch_idx = 'OS-{}'.format(sys.argv[2][-6:-4])

figname = '{}-OS-{}.png'.format(ccd_idx, ch_idx)

def f(a,x):
    """
    f = a[0]*x + a[1]
    """
    return a[0]*x + a[1]

prc_model = odr.Model(f)

tab = Table.read(sys.argv[2], format='ipac')
t = np.array(tab['t'])
flux = np.array(tab['dn_bin'])
tt = t - t.min()*0
fflux = flux - flux.min()*0

#print(t)

gain = 1.55

flux_max = np.max(fflux)
idx = np.logical_and( fflux > 0.02*flux_max, fflux < 0.8*flux_max)

fig = plt.figure(figsize=(8,6))

# fit prc
tfit, xfit = tt[idx], fflux[idx]
prc_data = odr.Data(tfit, xfit, wd=1./np.power(tfit,2), we=1/np.power(xfit,2))
#prc_data = odr.Data(tfit, xfit, wd=1., we=1)
#fitter = odr.ODR(prc_data, prc_model, ifixb=[1,0], beta0=[0,np.median(xfit/tfit)])
fitter = odr.ODR(prc_data, prc_model, ifixb=[0,0], beta0=[0,np.median(xfit/tfit)])
res = fitter.run()
#res.pprint()
k = res.beta[1]

# non-lin

nonlin = (fflux[1:] - tt[1:]*k) / (tt[1:]*k)

# plot PRC

ax1 = fig.add_subplot(2,1,1)
plt.plot(tt, tt*k/1000, 'm--', label='linear fitting')
plt.scatter(tt, fflux/1000, marker='o', c='y', alpha=0.75)
plt.scatter(tfit, xfit/1000, marker='o', c='g', label='used data:')
plt.xlim(-0.1,tt.max()*1.05)
plt.xlabel('Time(s)')
plt.ylabel('Flux (k ADU)')
plt.legend()

plt.title('{}, {}'.format(ccd_idx,ch_idx))

ax2 = fig.add_subplot(2,1,2)

plt.hlines(+0.018*100, xmin=0, xmax=1.1*tt.max(), linestyle='dashed', color='grey', label=r'$+1.8$')
plt.hlines(0, xmin=0, xmax=1.1*tt.max(), linestyle='dashed', color='grey')
plt.hlines(-0.018*100, xmin=0, xmax=1.1*tt.max(), linestyle='dashed', color='grey', label=r'$-1.8$')
plt.scatter(tt[1:], nonlin*100, marker='o')
plt.xlim(-0.1,tt.max()*1.05)
plt.xlabel('Time(s)')
plt.ylabel('Response Linearity (%)')
plt.legend()

plt.subplots_adjust(top=0.95, bottom=0.085, hspace=0.2)

plt.savefig(figname)

#plt.show()
